---
output: html_document
editor_options: 
  chunk_output_type: console
---

# Load packages and data

```{r packages and set wd}

# load packages
pacman::p_load(tidyverse, zoo, lubridate, tidyquant, reshape,
               pscl, cowplot, texreg, stargazer, broom)

`%nin%` = Negate(`%in%`)

theme <- theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position="top",
        strip.background = element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(),
        panel.grid.major= element_line(color = "seashell1"),
        plot.background=element_blank(),
        axis.text.x = element_text(size = 8),
        legend.text = element_text(size = 8),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_text(size = 8))

gov_reps = c("Mitchell_Grant","Bellemare_Diane","Harder_Peter")

```

```{r}

month_data <- readRDS("./data/transformed/month_data.RDS")
 
# Get rid of a few duplicates where people change party mid month
month_data = month_data %>% 
  ungroup() %>% 
  group_by(year_month, name_parls, senator, gender, province, treated, 
           appointing_party, appointing_PM, appointed_after, start, 
           og_isg, con_affiliation, lib_affiliation, independent, partisan,
           ind, year, trudeau_gov, ISG_lead, ISG_lag, ISG, ISG_linear, 
           IAB_lead, IAB_lag, IAB, IAB_linear, tenure_years, log_tenure, senior) %>% 
  summarize(party = case_when(
    str_detect(paste(party, collapse = "|"),"ISG") ~ "ISG",
    str_detect(paste(party, collapse = "|"),"NONE") ~ "NONE",
    str_detect(paste(party, collapse = "|"),"LPC") ~ "LPC",
    str_detect(paste(party, collapse = "|"),"CPC") ~ "CPC",
    str_detect(paste(party, collapse = "|"),"NDP") ~ "NDP",
    str_detect(paste(party, collapse = "|"),"BQ") ~ "BQ",
    str_detect(paste(party, collapse = "|"),"GREEN") ~ "Green"),
    total = sum(total),
    member = max(member),
    leader = max(leader),
    total_sens = max(total_sens))

month_data = month_data %>% 
  ungroup() %>% 
  mutate(gov_rep = factor(ifelse(as.character(name_parls) %in% gov_reps, 1, 0)))

num_independent = month_data %>%
  filter(senator == 1) %>% 
  ungroup() %>% 
  group_by(year_month) %>% 
  summarize(
    total_nonpartisan = sum(ifelse(party %in% c("NONE","ISG"),1,0)),
    num_ind = sum(independent, na.rm = TRUE),
    n = n()
  ) %>% 
  mutate(
    prop_ind = total_nonpartisan/n
  )

month_data = merge(month_data, num_independent, by = "year_month")

month_data = month_data %>% 
  mutate(interim = prop_ind - min(prop_ind),
         prop_ind_rescaled = interim/max(interim))

month_data <- month_data %>%
  mutate(period = (year_month - min(month_data$year_month))*12 + 1,
         IAB_lead = ifelse(year_month >= as.yearmon("Oct 2015") & 
                             year_month < as.yearmon("Jan 2016"), 1, 0),
         ISG_lead = ifelse(year_month >= as.yearmon("Dec 2015") & 
                             year_month < as.yearmon("Mar 2016"), 1, 0))

month_data <- month_data %>%
  mutate(mp_treatment = case_when(
           lib_affiliation == 1 & senator == 0 ~ 1,
           lib_affiliation == 0 & senator == 0 ~ -1,
           senator == 1 ~ 0
         ))

month_data <- month_data %>%
  mutate(first_year = case_when(
    year_month >= as.yearmon("Oct 2015") & year_month < as.yearmon("Nov 2016") ~ 1,
    TRUE ~ 0)
  )

month_data = month_data %>% 
  mutate(category = factor(case_when(
    senator == 1 & independent == 1 ~ "ind",
    senator == 1 & independent == 0 ~ "sen",
    senator == 0 ~ "MP",
  ), levels = c("MP","sen","ind")))

month_data = mutate(month_data, posts = leader + member)

```

```{r}

month_data %>%
  group_by(year_month) %>%
  summarize(mean = sum(total)) %>% 
  ggplot(aes(x = year_month, y= mean)) +
  geom_bar(stat = "identity") + theme_bw()

```

```{r}

month_data %>% 
  filter(senator == 1) %>% 
  group_by(gov_rep, year_month) %>% 
  summarize(total = sum(total)/n()) %>% 
  ggplot(aes(x = year_month, y = total)) +
  geom_bar(stat = "identity") +
  facet_wrap(gov_rep~.)

```

# Figure 1

```{r}

fig1_dat <- filter(month_data, !is.na(senator)) %>%
  filter(year_month >= "Aug 2010") %>%
  group_by(year_month, senator) %>%
  dplyr::summarise(total = sum(total),
                   n = n()) %>%
  mutate(senator = factor(senator,
                          labels = c("Members of Parliament","Senators"),
                          levels = c(0,1)),
         total = total/n)

fig1 <- fig1_dat %>%
  ggplot(aes(x = year_month, y = total, col = senator, fill = senator)) +
  geom_col(position = "stack") + 
  geom_ma(size = 2, n = 7, linetype = "solid") + 
  scale_x_yearmon(format = "%b %Y", n = 12) +
  scale_fill_manual(name = "",
                    breaks = c("Members of Parliament","Senators"),
                    values = c("white","white")) +
  scale_color_manual(name = "",
                     breaks = c("Members of Parliament","Senators"),
                     values = c("slategrey","black")) +
  theme +
  labs(x = "", y = "Average number of contacts between\nparliamentarians and lobbyists per month") +
  facet_grid(senator~.) +
  geom_vline(xintercept = as.yearmon("Jan 2016"), linetype = "dashed") +
  guides(fill = FALSE, color = FALSE)

fig1_text <- data.frame(senator = factor("Members of Parliament", 
                                        levels = c("Members of Parliament",
                                                   "Senators")))
fig1 + geom_text(data = fig1_text,
              mapping = aes(x = as.yearmon("Feb 2015"),
                    y = 3.5, 
                    label = "IAB established"),
              size = 3,
              color = "black",
              show.legend=FALSE) +
  coord_cartesian(ylim=c(0,4)) + 
  theme(legend.position="top")

ggsave("outputs/figures/fig1.png", height = 5, width = 6.5)
ggsave("outputs/figures/fig1.eps", height = 5, width = 6.5)

```

# Table 1

```{r}

tab1 = list()

tab1[[1]] = pscl::zeroinfl(data = month_data, dist = "negbin",
                 total ~ tenure_years + gender + member*senator + leader*senator + IAB*senator)

tab1[[2]] = pscl::zeroinfl(data = month_data, dist = "negbin",
                 total ~ tenure_years + gender + member*senator + leader*senator + prop_ind_rescaled*senator)

tab1[[3]] = pscl::zeroinfl(data = month_data, dist = "negbin",
                 total ~ tenure_years + gender + member + leader + IAB*category)

tab1[[4]] = pscl::zeroinfl(data = month_data, dist = "negbin",
                 total ~ tenure_years + gender + member + leader + prop_ind_rescaled*category)

screenreg(tab1, single.row = TRUE, omit.coef = "(year_month)|(name_parl)|(province)")

texreg(tab1,
       file = "outputs/tables/tab1.tex", dcolumn = TRUE, use.packages = FALSE,
       scalebox = 0.65, float.pos = "ht!", label = "table:tab1", single.row = TRUE,
       custom.model.names = c("1: H1A","2: H1B","3: H2A","4: H2B"),
       omit.coef = "(year_month)|(name_parl)|(province)|(theta)|(tenure_years)|(gender)|(member)|(leader)|(Intercept)|(period)",
       caption = "Lobbying activity directed at Canadian Parliamentarians", caption.above = TRUE,
       custom.note = ("\\parbox{1.5\\linewidth}{\\vspace{2pt}%stars. Zero-inflated Negative Binomial Regression models showing log-odds with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month.}"),
          custom.coef.names = c("Senator",
                                "Post-IAB",
                                "Post-IAB x Senator",
                                "Z: Senator",
                                "Z: Post-IAB",
                                "Z: Post-IAB x Senator",

                                "Proportion Independent",
                                "Senator x Proportion Independent",
                                "Z: Proportion Independent",
                                "Z: Senator x Proportion Independent",

                                "Partisan Senator",
                                "Independent Senator",
                                "Post-IAB x Partisan Senator",
                                "Post-IAB x Independent Senator",
                                "Z: Partisan Senator",
                                "Z: Independent Senator",
                                "Z: Post-IAB x Partisan Senator",
                                "Z: Post-IAB x Independent Senator",

                                "Proportion Independent x Partisan Senator",
                                "Proportion Independent x Independent Senator",
                                "Z: Proportion Independent x Partisan Senator",
                                "Z: Proportion Independent x Independent Senator"),
       groups = list("DID House versus Senate" = 1:6,
                     "Proportion House versus Senate" = 7:10,
                     "DID Parliamentary role (relative to MPs)" = 11:18,
                     "Proportion Parliamentary role (relative to MPs)" = 19:22),
       reorder.coef = c(1:6, 7:10, 11:18, 19:22))

```

# Figure 2

## Left panel

```{r}

predict_dat = filter(month_data) %>% 
  select(tenure_years, gender, member, leader, category, IAB) %>% 
  mutate(id = row_number())

predict_dat = rbind(predict_dat %>% mutate(instance = 1),
                    predict_dat %>% mutate(instance = 2),
                    predict_dat %>% mutate(instance = 3),
                    predict_dat %>% mutate(instance = 4),
                    predict_dat %>% mutate(instance = 5),
                    predict_dat %>% mutate(instance = 6))

predict_dat = mutate(predict_dat,
                     IAB = ifelse(instance %in% c(1,2,3), 1, 0),
                     category = case_when(
                       instance %in% c(1,4) ~ "MP",
                       instance %in% c(2,5) ~ "ind",
                       instance %in% c(3,6) ~ "sen"
                     ))

predict_dat[,c("fit")] <- predict(tab1[[3]], type = "response", newdata = predict_dat)

predict_dat %>% 
  group_by(IAB, category) %>% 
  summarize(mean = mean(fit), n = n())

ggdat = predict_dat %>% 
  select(instance, id, IAB, category, fit) %>% 
  mutate(characteristics = case_when(
    instance == 1 ~ "post_mp",
    instance == 2 ~ "post_ind",
    instance == 3 ~ "post_sen",
    instance == 4 ~ "pre_mp",
    instance == 5 ~ "pre_ind",
    instance == 6 ~ "pre_sen"
  )) %>% 
  select(characteristics, id, fit) %>% 
  pivot_wider(id_cols = id, names_from = characteristics, values_from = fit) %>% 
  mutate(mps = (post_mp/pre_mp)-1, sens = (post_sen/pre_sen)-1, inds = (post_ind/pre_ind) - 1) %>% 
  select(id, mps, sens, inds) %>% 
  gather(key = "key", value = "value", - id) %>% 
  group_by(key) %>% 
  summarize(mean = mean(value), sd = sd(value), n = n()) %>% 
  mutate(se = sd/sqrt(n)) %>% 
  mutate(key = ordered(
    case_when(
      key == "inds" ~ "Independent\nSenators",
      key == "sens" ~ "Partisan\nSenators",
      key == "mps" ~ "Members of\nParliament"),
    levels = c("Members of\nParliament","Partisan\nSenators","Independent\nSenators")))

theme_update(legend.key = element_rect(size = 6, fill = "white", colour = NA),
             legend.key.size = unit(1, "cm"))

p1 = ggdat %>% 
  ggplot(aes(x = key, y = mean)) +
  geom_point(aes(colour = key, shape = key, fill = key), size = 5, color = "black") +
  theme +
  scale_fill_brewer(type = "div", palette = "Greys") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1.5), breaks = seq(0,1.5,0.25)) +
    scale_shape_manual(values = c(21,22,23)) + 
  labs(x = "", y = "Average increase in contacts per month after IAB") +
  theme(legend.title = element_blank()) +
  guides(col = guide_legend(ncol = 3,keyheight = 2)); p1
```

## Right panel

```{r}

right_model = pscl::zeroinfl(data = month_data, dist = "negbin",
                 total ~ tenure_years + gender + member + leader + prop_ind*category)


predict_dat = filter(month_data) %>% 
  select(tenure_years, gender, member, leader, category, IAB) %>% 
  mutate(category = factor(category)) %>% 
  mutate(id = row_number())

prop_ind_range = seq(min(month_data$prop_ind), max(month_data$prop_ind), 0.01)

predict_dat2 = expand.grid.df(predict_dat,
                              data.frame(prop_ind = prop_ind_range,
                                         instance = seq(1, length(prop_ind_range))))

predict_dat2 = rbind(predict_dat2 %>% 
                       mutate(category = unique(predict_dat$category)[1]),
                     predict_dat2 %>% 
                       mutate(category = unique(predict_dat$category)[2]),
                     predict_dat2 %>% 
                       mutate(category = unique(predict_dat$category)[3]))

predict_dat2[,c("fit")] <- predict(right_model, type = "response",
                                   newdata = predict_dat2)

ggdat2 = predict_dat2 %>% 
  group_by(category, prop_ind) %>% 
  summarize(mean = mean(fit), sd = sd(fit), n = n()) %>% 
  mutate(se = sd/sqrt(n)) %>% 
  ungroup() %>% 
  mutate(key = ordered(
    case_when(
      category == "ind" ~ "Independent\nSenators",
      category == "sen" ~ "Partisan\nSenators",
      category == "MP" ~ "Members of\nParliament"),
    levels = c("Members of\nParliament","Partisan\nSenators","Independent\nSenators")))

p2 = ggdat2 %>% 
  ggplot(aes(x = prop_ind, y = mean, color = key, fill = key, linetype = key)) +
  geom_point(aes(colour = key, shape = key), size = 2, col = "black") +
  theme + 
  scale_y_continuous(limits = c(0, 1.6)) + 
  scale_fill_brewer(type = "div", palette = "Greys") +
  scale_shape_manual(values = c(21,22,23)) + 
  guides(color = FALSE, shape = FALSE, fill = FALSE) + 
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(x = "Percentage of Senators that are independent", y = "Average number of contacts per month"); p2

```

## Save


```{r}

prow <- plot_grid(
  p1 + theme(legend.position="none"),
  p2 + theme(legend.position="none"),
  align = 'vh',
  hjust = -1,
  nrow = 1
)

legend <- get_legend(
  p1 + 
    theme(legend.position = "top")
)

# add the legend to the row we made earlier. Give it one-third of 
# the width of one plot (via rel_widths).
plot_grid(legend, prow, ncol = 1, rel_widths = c(3, 0.1), rel_heights = c(0.3,3))

ggsave("outputs/figures/fig2.eps", height = 5, width = 6.5)

```

# Appendix

## A: Control Variables

```{r}

texreg(tab1,
       file = "outputs/tables/tab1_app.tex", dcolumn = TRUE, use.packages = FALSE,
       scalebox = 0.65, float.pos = "ht!", label = "table:tab1_app", single.row = TRUE,
       custom.model.names = c("1: H1A","2: H1B","3: H2A","4: H2B"),
       omit.coef = "(year_month)|(name_parl)|(province)|(theta)",
       caption = "Lobbying activity directed at Canadian Parliamentarians", caption.above = TRUE,
       custom.note = ("\\parbox{1.5\\linewidth}{\\vspace{2pt}%stars. Zero-inflated Negative Binomial Regression models showing log-odds with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month. All models include province fixed effects.}"),
          custom.coef.names = c("Intercept",
                                "Tenure (Years)",
                                "Male",
                                "No. of committee memberships",
                                "Senator",
                                "No. of committee leaderships",
                                "Post-IAB",
                                "Senator x No. of committee memberships",
                                "Senator x No. of committee leaderships",
                                "Post-IAB x Senator",

                                "Z: Intercept",
                                "Z: Tenure (Years)",
                                "Z: Male",
                                "Z: No. of committee memberships",
                                "Z: Senator",
                                "Z: No. of committee leaderships",
                                "Z: Post-IAB",
                                "Z: Senator x No. of committee memberships",
                                "Z: Senator x No. of committee leaderships",
                                "Z: Post-IAB x Senator",

                                "Proportion Independent",
                                "Senator x Proportion Independent",
                                "Z: Proportion Independent",
                                "Z: Senator x Proportion Independent",
                                
                                "Partisan Senator",
                                "Independent Senator",
                                "Post-IAB x Partisan Senator",
                                "Post-IAB x Independent Senator",
                                "Z: Partisan Senator",
                                "Z: Independent Senator",
                                "Z: Post-IAB x Partisan Senator",
                                "Z: Post-IAB x Independent Senator",
                                
                                "Proportion Independent x Partisan Senator",
                                "Proportion Independent x Independent Senator",
                                
                                "Z: Proportion Independent x Partisan Senator",
                                "Z: Proportion Independent x Independent Senator"))

```

## B: Descriptive statistics

```{r}

desc1 <- filter(month_data) %>% 
  ungroup() %>%
  mutate(gender = ifelse(gender == "M", 1, 0),
         year = as.numeric(as.character(year))) %>%
  dplyr::select(total, tenure_years, gender, 
                member, leader, senator, prop_ind) %>%
  psych::describe(skew = TRUE, check = TRUE) %>%
  dplyr::select(n, mean, sd, median, min, max) %>%
  rownames_to_column()

colnames(desc1) <- c("Variable","n","Mean","SD","Median","Min","Max")

desc1[,'Variable'] <- c("Number of contacts in a month",
                        "Tenure (Years)",
                        "Male",
                        "Number of committee memberships",
                        "Number of committee leaderships",
                        "Senator",
                        "Proportion Independent")

print(xtable::xtable(desc1, 
                     digits = c(0,0,0,2,2,2,0,0),
                     label = "table:descriptives_hoc",
                     caption = "Descriptives for full parliamentarian sample"),
      caption.placement = "top",
      file="outputs/tables/descriptives_hoc.tex")

```

```{r}

desc2 <- filter(month_data, senator == 1) %>% 
  ungroup() %>%
  mutate(gender = ifelse(gender == "M", 1, 0),
         year = as.numeric(as.character(year))) %>%
  dplyr::select(total, tenure_years, gender, 
                member, leader, independent, prop_ind) %>%
  psych::describe(skew = TRUE, check = TRUE) %>%
  dplyr::select(n, mean, sd, median, min, max) %>%
  rownames_to_column()

colnames(desc2) <- c("Variable","n","Mean","SD","Median","Min","Max")
desc2[,'Variable'] <- c("Number of contacts in a month",
                        "Tenure (Years)",
                        "Male",
                        "Number of committee memberships",
                        "Number of committee leaderships",
                        "Independent",
                        "Proportion Independent")

print(xtable::xtable(desc2,
                     digits = c(0,0,0,2,2,2,0,0),
                     label = "table:descriptives_sen",
                     caption = "Descriptives for Senate sample"), 
      caption.placement = "top",
      file="outputs/tables/descriptives_sen.tex")

```

## C: DID

### Parallel trends

```{r}

parallel_trend <- list()

parallel_trend[[1]] <- pscl::zeroinfl(data = filter(month_data, IAB == 0), dist = "negbin",
                                 total ~ tenure_years + gender + province + 
                                   member*senator + leader*senator +
                                   period*senator)

parallel_trend[[2]] <- pscl::zeroinfl(data = filter(month_data, IAB == 0), dist = "negbin",
                                 total ~ tenure_years + gender + province + 
                                   member + leader + period*category)

screenreg(parallel_trend, single.row = TRUE, omit.coef = "(year_month)|(name_parl)|(province)")

# Very slight
texreg(parallel_trend, 
          scalebox = 0.8,
          float.pos = "ht!",
          label = "table:parallel",
          file = "outputs/tables/parallel.tex",
          caption = "Testing parallel pre-trend of lobbying activity directed at Canadian Parliamentarians",
          caption.above = TRUE,
          dcolumn = TRUE,use.packages = FALSE,
          single.row = TRUE,
          custom.note = ("\\parbox{1.1\\linewidth}{\\vspace{2pt}%stars. Zero-inflated negative binomial regression models showing log-odds with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month. All models on subset of data that excludes all observations after the establishment of the IAB.}"),
          omit.coef = "(name_parl)|(province)|(theta)",
          custom.coef.names = c("Intercept",
                                "Tenure (Years)",
                                "Male",
                                "No. of committee memberships",
                                "Senator",
                                "No. of committee leaderships",
                                "Trend",
                                "No. of committee memberships x Senator",
                                "No. of committee leaderships x Senator",
                                "Trend x Senator",
                                
                                "Z: Intercept",
                                "Z: Tenure (Years)",
                                "Z: Male",
                                "Z: No. of committee memberships",
                                "Z: Senator",
                                "Z: No. of committee leaderships",
                                "Z: Trend",
                                "Z: No. of committee memberships x Senator",
                                "Z: No. of committee leaderships x Senator",
                                "Z: Trend x Senator",

                                "Partisan Senator",
                                "Independent Senator",
                                "Trend x Partisan Senator",
                                "Trend x Independent Senator",
                                
                                "Z: Partisan Senator",
                                "Z: Independent Senator",
                                "Z: Trend x Partisan Senator",
                                "Z: Trend x Independent Senator"),
          custom.model.names = c("1: HoC versus Senator","2: Parliamentary Groups"))

```

### Leads

```{r}

leads <- list()

leads[[1]] <- pscl::zeroinfl(data = filter(month_data), dist = "negbin",
                                 total ~ tenure_years + gender + province + 
                                   member*senator + leader*senator + IAB*senator + IAB_lead*senator)

leads[[2]] <- pscl::zeroinfl(data = filter(month_data), dist = "negbin",
                                 total ~ tenure_years + gender + province + 
                                   member + leader + IAB*category + IAB_lead*category)

screenreg(leads, omit.coef = "(year_month)|(name_parl)|(province)", single.row = TRUE)

texreg(leads, 
       scalebox = 0.8,
       float.pos = "ht!",
       label = "table:leads",
       file = "outputs/tables/leads.tex",
       caption = "Testing leads for lobbying activity directed at Canadian Parliamentarians",
       caption.above = TRUE,
       dcolumn = TRUE,use.packages = FALSE,
       single.row = TRUE,
       custom.note = ("\\parbox{1.1\\linewidth}{\\vspace{2pt}%stars. Zero-inflated Negative Binomial Regression models showing log-odds with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month. All models include province fixed effects.}"),
       omit.coef = "(name_parl)|(province)|(theta)",
       custom.coef.names = c("Intercept",
                             "Tenure (Years)",
                             "Male",
                             "No. of committee memberships",
                             "Senator",
                             "No. of committee leaderships",
                             "Post-IAB",
                             "IAB Lead",
                             "No. of committee memberships x Senator",
                             "No. of committee leaderships x Senator",
                             "Post-IAB x Senator",
                             "IAB Lead x Senator",
                             
                             "Z: Intercept",
                             "Z: Tenure (Years)",
                             "Z: Male",
                             "Z: No. of committee memberships",
                             "Z: Senator",
                             "Z: No. of committee leaderships",
                             "Z: Post-IAB",
                             "Z: IAB Lead",
                             "Z: No. of committee memberships x Senator",
                             "Z: No. of committee leaderships x Senator",
                             "Z: Post-IAB x Senator",
                             "Z: IAB Lead x Senator",
                             
                             "Partisan Senator",
                             "Independent Senator",
                             "Post-IAB x Partisan Senator",
                             "Post-IAB x Independent Senator",
                             "IAB Lead x Partisan Senator",
                             "IAB Lead x Independent Senator",
                             
                             
                             "Z: Partisan Senator",
                             "Z; Independent Senator",
                             "Z; Post-IAB x Partisan Senator",
                             "Z: Post-IAB x Independent Senator",
                             "Z: IAB Lead x Partisan Senator",
                             "Z: IAB Lead x Independent Senator"),
       custom.model.names = c("1: HoC versus Senator","2: Parliamentary Groups"))
       
```

### Exclusion

#### Visualization

```{r}

filter(month_data, senator == 1) %>% 
  group_by(year_month, party, independent) %>%
  summarize(n = n()) %>% 
  mutate(independents = ifelse(party %in% c("NONE","ISG"), 1, 0 )) %>%
  mutate(treated = factor(independent,
                          ordered = TRUE,
                          labels = c("Party-affiliated","Independent"),
                          levels = c(0,1))) %>%
  arrange(desc(n)) %>% 
  ggplot(aes(x = year_month, y = n, fill = treated)) +
  geom_col(position = "stack") + 
  theme + 
  scale_x_yearmon(format = "%b %Y", n = 12) +
  scale_fill_manual(name = "",
                    breaks = c("Party-affiliated","Independent"),
                    values = c("lightgrey","black")) + 
  labs(x = "", y = "Number of Senators")

ggsave("outputs/figures/independents.png", height = 5, width = 6.5)

```

### Models

```{r}

excl <- list()

excl[[1]] <- pscl::zeroinfl(data = filter(month_data, senator == 1, og_isg == 0), dist = "negbin",
                            total ~ tenure_years + gender + member + leader + independent*IAB)

excl[[2]] <- pscl::zeroinfl(data = filter(month_data, senator == 1, independent == 0 | og_isg == 1), dist = "negbin",
                            total ~ tenure_years + gender + member + leader + independent*IAB)

excl[[3]] <- pscl::zeroinfl(data = filter(month_data, senator == 1, appointed_after == 0), dist = "negbin",
                            total ~ tenure_years + gender + member + leader + independent*IAB)

screenreg(excl, single.row = TRUE, omit.coef = "(year_month)|(name_parl)|(province)")

texreg(excl,
       scalebox = 0.8,
       float.pos = "ht!",
       label = "table:sen_exclusion",
       file = "outputs/tables/sen_exclusion.tex",
       caption = "Lobbying activity directed at Canadian Senators (subsamples)",
       caption.above = TRUE,
       dcolumn = TRUE,
       use.packages = FALSE,
       single.row = TRUE,
       custom.note = ("\\parbox{1.1\\linewidth}{\\vspace{2pt}%stars. Zero-inflated Poisson Regressions with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month. All models run on subsamples of Senator data and include province fixed effects.}"),
       custom.model.names = c("1: No orig. ISG","2: Orig. ISG",
                              "3: Pre-IAB App."),
       
       custom.coef.names = c("Intercept","Tenure (Years)","Male",
                             "No. of committee memberships",
                             "No. of committee leaderships",
                             "Independent",
                             "Post-IAB",
                             "Independent x Post-IAB",
                             
                             "Z: Intercept",
                             "Z: Tenure (Years)",
                             "Z: Male",
                             "Z: No. of committee memberships",
                             "Z: No. of committee leaderships",
                             "Z: Independent",
                             "Z: Post-IAB",
                             "Z: Independent Joiner x Post-IAB"),
       groups = list("Count model" = 1:8, "Zero model" = 9:16),
       omit.coef = "(year_month)|(name_parl)|(province)|(theta)")

```

## Change in Senators

```{r}

month_data_alt <- month_data %>% 
  mutate(independent = ifelse(name_parls %in% gov_reps, 0, independent))

num_independent_alt = month_data_alt %>%
  filter(senator == 1) %>% 
  ungroup() %>% 
  group_by(year_month) %>% 
  summarize(
    total_nonpartisan = sum(ifelse(party %in% c("NONE","ISG") & gov_rep != 1, 1, 0)),
    num_ind = sum(independent, na.rm = TRUE),
    n = n()
  ) %>% 
  mutate(
    prop_ind = total_nonpartisan/n
  )

month_data_alt = merge(month_data_alt, num_independent_alt, by = "year_month")

month_data_alt = month_data_alt %>% 
  mutate(category = factor(case_when(
    senator == 1 & independent == 1 ~ "ind",
    senator == 1 & independent == 0 ~ "sen",
    senator == 0 ~ "MP",
  ), levels = c("MP","sen","ind")))

gov_reps = list()

gov_reps[[1]] = pscl::zeroinfl(data = month_data_alt, dist = "negbin",
                 total ~ tenure_years + gender + member*senator + leader*senator + IAB*senator)

gov_reps[[2]] = pscl::zeroinfl(data = month_data_alt, dist = "negbin",
                 total ~ tenure_years + gender + member*senator + leader*senator + prop_ind_rescaled*senator)

gov_reps[[3]] = pscl::zeroinfl(data = month_data_alt, dist = "negbin",
                 total ~ tenure_years + gender + member + leader + IAB*category)

gov_reps[[4]] = pscl::zeroinfl(data = month_data_alt, dist = "negbin",
                 total ~ tenure_years + gender + member + leader + prop_ind_rescaled*category)

screenreg(gov_reps, single.row = TRUE,
          omit.coef = "(year_month)|(name_parl)|(province)|(theta)|(tenure_years)|(gender)|(member)|(leader)|(Intercept)|(period)")

texreg(gov_reps,
       file = "outputs/tables/gov_reps.tex", dcolumn = TRUE, use.packages = FALSE,
       scalebox = 0.65, float.pos = "ht!", label = "table:gov_reps", single.row = TRUE,
       custom.model.names = c("1: H1A","2: H1B","3: H2A","4: H2B"),
       omit.coef = "(year_month)|(name_parl)|(province)|(theta)|(tenure_years)|(gender)|(member)|(leader)|(Intercept)|(period)",
       caption = "Lobbying activity directed at Canadian Parliamentarians", caption.above = TRUE,
       custom.note = ("\\parbox{1.5\\linewidth}{\\vspace{2pt}%stars. Zero-inflated Negative Binomial Regression models showing log-odds with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month. All models include province fixed effects.}"),
          custom.coef.names = c("Senator",
                                "Post-IAB",
                                "Post-IAB x Senator",
                                "Z: Senator",
                                "Z: Post-IAB",
                                "Z: Post-IAB x Senator",

                                "Proportion Independent",
                                "Senator x Proportion Independent",
                                "Z: Proportion Independent",
                                "Z: Senator x Proportion Independent",

                                "Partisan Senator",
                                "Independent Senator",
                                "Post-IAB x Partisan Senator",
                                "Post-IAB x Independent Senator",
                                "Z: Partisan Senator",
                                "Z: Independent Senator",
                                "Z: Post-IAB x Partisan Senator",
                                "Z: Post-IAB x Independent Senator",

                                "Proportion Independent x Partisan Senator",
                                "Proportion Independent x Independent Senator",
                                "Z: Proportion Independent x Partisan Senator",
                                "Z: Proportion Independent x Independent Senator"),
       groups = list("DID House versus Senate" = 1:6,
                     "Proportion House versus Senate" = 7:10,
                     "DID Parliamentary role (relative to MPs)" = 11:18,
                     "Proportion Parliamentary role (relative to MPs)" = 19:22),
       reorder.coef = c(1:6, 7:10, 11:18, 19:22))

```

## Linear models

### Fixed effects

```{r}


lms2 <- list()

lms2[[1]] <- lm(data = month_data, 
                 total ~ tenure_years + gender + province +
                   member*senator + leader*senator + 
                 factor(year_month) + factor(name_parls) + IAB*senator)

lms2[[2]] <- lm(data = month_data, 
                 total ~ tenure_years + gender + province +
                   member*senator + leader*senator + 
                 factor(year_month) + factor(name_parls) + senator*prop_ind_rescaled)

lms2[[3]] <- lm(data = month_data, 
                 total ~ province + tenure_years + gender + province +
                 factor(year_month) + factor(name_parls) +
                   member + leader + IAB*category)

lms2[[4]] <- lm(data = month_data, 
                 total ~ province + tenure_years + gender + province +
                 factor(year_month) + factor(name_parls) + 
                   member + leader + category*prop_ind_rescaled)

screenreg(lms2, omit.coef = "(year_month)|(name_parl)|(province)",single.row = TRUE)

texreg(lms2, 
          #single.row = TRUE, omit.coef = "(year_month)|(name_parl)|(province)")
          scalebox = 0.6,
          float.pos = "ht!",
          label = "table:lms2",
          caption.above = TRUE,
          dcolumn = TRUE,
       use.packages = FALSE,
          single.row = TRUE,
          file = "outputs/tables/lms2.tex",
          caption = "Linear models of main-paper models (with individual and period fixed effects)",
          custom.note = ("\\parbox{1.3\\linewidth}{\\vspace{2pt}%stars. Linear models with standard errors in parentheses. Dependent variable: number of contacts with lobbyists per parliamentarian month. All models include province, period, and individual fixed-effects.}"),
          omit.coef = "(year_month)|(name_parl)|(province)",
       custom.coef.names = c("Intercept",
                             "Tenure (Years)",
                             "Male",
                             "No. of committee memberships",
                             "Senator",
                             "No. of committee leaderships",
                             "Senator x No. of committee memberships",
                             "Senator x No. of committee leaderships",
                             "Senator x Post-IAB",
                             "Senator x Proportion independent",
                             
                             "Partisan Senator x Post-IAB",
                             "Independent Senator x Post-IAB",
                             "Partisan Senator x Independent proportion",
                             "Independent Senator x Independent proportion"),
          custom.model.names = c("1: HoC vs Senate DID",
                                 "2: HoC vs Senate Proportion",
                                 "3: Grouping DID",
                                 "4: Grouping Proportion"))
```

```{r}

month_data %>% 
  group_by(name_parls, appointing_PM) %>% 
  count() %>% 
  group_by(appointing_PM) %>% 
  count()
```

